Journal of Theoretical Biology

Neutrophil swarming is a critical immune response in mammals and fish, in which neutrophils are recruited to inflammatory sites where they coordinate into a swarm that neutralizes pathogens. While excessive swarming can drive prolonged inflammation, a quantitative understanding of swarming dynamics remains limited. We developed a one-dimensional radial reaction-diffusion model of neutrophil swarming with two kinetic parameters, in order to capture the self-limiting swarming dynamics in both murine and human neutrophils in response to different inflammatory stimulus sizes. To ensure that the inverse problem is well-posed, we first performed sensitivity and identifiability analyses. We then developed a physics-informed neural network (PINN) to infer the key parameters governing swarm expansion and self-limitation. To account for uncertainty in noisy experimental measurements, we further extended this framework to a Bayesian PINN (B-PINN), which provides credible intervals for the inferred parameters. Both models were validated against synthetic data generated by numerical simulation and subsequently applied to in vitro experimental data from human and murine neutrophils in response to three bioparticle cluster sizes. The PINN-inferred dynamics show that larger bioparticle clusters are associated with greater cumulative recruitment and larger swarms in both species. The models further reveal species-specific differences in both the amplitude of initial recruitment and the timescale on which it self-limits. Additionally, the B-PINN posterior distributions quantify uncertainty in these species- and cluster size-dependent trends and identify where additional measurements would be most informative. To our knowledge, this is the first application of physics-informed machine learning to model neutrophil swarming dynamics. This framework provides a starting point for systematically comparing recruitment dynamics between human and murine neutrophils and offers guidance for future experimental design.

Biological populations have been evolving in multi-dimensional trait spaces. Ecological interactions can induce the evolutionary diversification of a biological population into distinct multiple populations, which is called evolutionary branching in adaptive dynamics theory. If possible mutations are restricted to certain directions due to developmental, physiological, or physical constraints, the resulting evolutionary dynamics may be restricted to subspaces with fewer dimensions than the original trait spaces. In general, such constraints are also evolving slowly, due to rare and slight mutational deviations from the constraint subspaces. In this sense, constraints are not strict. However, the method for describing non-strict constraints that is globally consistent has yet to be developed. This study proposes a general and flexible method for constructing a non-strict constraint subspace from the corresponding strict constraint subspace. This method enables efficient analysis of adaptive evolution including evolutionary branching under non-strict constraints in the framework of adaptive dynamics theory. Applying this method to a resource competition model shows that the shape of the carrying capacity distribution can be treated as a non-strict constraint curve, and that the likelihood of evolutionary branching is affected not only by the shape of the constraint curve but also by its slow evolution.

Ecological interactions play a crucial role in driving speciation and adaptive radiation. According to the adaptive dynamics theory of evolution in one-dimensional trait spaces, such evolutionary diversifications can occur when evolutionary branching points (convergence stable and evolutionarily unstable singular points) exist in the trait spaces. However, these evolutionary diversifications may occur across multiple traits at once, as exemplified by various within-species ecotypic diversifications involving multiple traits or genes. This paper attempt to extend this branching point condition to evolution in multi-dimensional trait spaces. Two criteria, termed branching possibility and branching inevitability, are developed to evaluate the likelihood of evolutionary branching in multi-dimensional trait spaces. Branching possibility guarantees the existence of evolutionary pathways that lead to branching (where deterministic parts of these pathways are described with the canonical equation of adaptive dynamics theory). Branching inevitability guarantees a steady progression of evolutionary branching (where the progression degree is described with a time-increasing function related to the local Lyapunov function, the potential function in an evolutionary potential game, and Fisher's fundamental theorem). Numerically simulated evolution in multi-dimensional trait spaces show that points having branching possibility or inevitability can induce evolutionary branching not only under asexual reproduction but also under sexual reproduction.

Tumor progression is often modeled as a passive response to external therapy or immune pressure, but tumor populations may also exhibit population-level regulation of proliferation and apoptosis. We develop a continuous-time Markov decision framework in which a controlled birth--death process represents a tumor population modulating the balance between proliferation and susceptibility to apoptosis in the presence of extrinsic death pressure. We examine threshold and quadratic costs, an unbounded linear reward, and constrained linear and quadratic formulations to determine how objective structure shapes optimal policies and induced population drift. Threshold and quadratic penalties generate restoring dynamics, with transitions from growth to suppression and regions of near-neutral drift associated with regulated or near-dormant behavior. An unbounded linear reward instead produces sustained or near-neutral growth without a restoring regime. Under constraints, a linear reward expands the region of positive drift as capacity increases, whereas a quadratic reward generates restoring, logistic-like drift around an interior population scale. These results show that regulated tumor dynamics depend on how growth incentives, extrinsic death pressure, penalties, and constraints scale with population size.

Despite the impacts of within-host disease dynamics on disease outcomes in individual hosts and disease spread among-hosts, generic models of within-host population dynamics have received far less attention than their among-host counterparts. While a number of models have been proposed to explore theoretical eco-evolutionary dynamics, they have yet to be evaluated for estimability, raising questions about their ability to provide reliable inference when confronted with data. We evaluated the estimability of two generic within-host population dynamics models by assessing: (1) parameter estimation, our ability to recover correct values of model parameters from data, (2) the consequences of mis-assigning the underlying mechanistic model on parameter estimation, and (3) the reproduction of qualitative dynamics, or, our ability to use parameter estimates to reproduce observed dynamical behaviours. In some cases, fitting a mis-matched mechanistic model to time series data produced reasonable parameter estimates that were able to reproduce system dynamics, and that when provided the data-generating model, parameter uncertainty can produce substantial behavioural uncertainty. Our findings highlight the impacts of structural, parametric, and behavioural uncertainty on inference, and demonstrate the value of improving system-specific knowledge to prevent the use of incorrect functional forms and of measuring consequential parameters to improve estimability.

Moment-closure methods provide a widely used approach for approximating low-order statistics of stochastic biochemical reaction networks, but the resulting closed equations may exhibit physically meaningless dynamics. Here, we systematically investigate the stability of second-order moment-closure approximations for a bursty gene expression model with nonlinear protein degradation. We consider the Gaussian, log-normal, and negative-binomial closures. By analyzing the resulting dynamical systems, we show that the Gaussian closure fails the proposed stability criterion for all parameter values, whereas both the log-normal and negative-binomial closures possess a unique physically accessible fixed point and are stable if and only if kB(1 + 2B) [≥] 1, where B is the mean burst size and k is the dimensionless protein synthesis rate. Moreover, this condition is automatically satisfied when the steady-state mean protein copy number is at least one, suggesting that it is mild in biologically relevant regimes. Numerical comparisons over a broad parameter range further show that the negative-binomial closure is the most accurate, while the Gaussian closure is the least accurate. We further develop higher-order negative-binomial moment-closure schemes and demonstrate that increasing the closure order substantially improves predictive accuracy. These results highlight the importance of considering dynamical stability, in addition to numerical accuracy, when assessing moment-closure approximations.

Fractional differential equations are increasingly used to represent temporal memory in tumour immune and pharmacokinetic models, but changing the derivative order without a consistent time scale can make comparisons difficult to interpret. We study a dimensionless, commensurate Caputo system for tumour, immune, normal cell, and drug states. We establish positivity, boundedness, the tumour invasion threshold, and local stability conditions. We then use a fractional Adams Bashforth Moulton predictor corrector method with an explicit mesh convergence study. Two conclusions correct potentially misleading interpretations of the model. First, the equilibrium locations are independent of the fractional order: decreasing the order changes transient time scales but does not by itself alter the model's asymptotic states. Second, after a unit bolus the drug equation has the exact solution C({tau})=E(-K{tau}). Its maximum is therefore the initial concentration, not a delayed peak. For 0<<1, the algebraic tail makes the infinite-horizon area under the concentration curve diverge. In the illustrative coupled simulations, the apparent late tumour reduction at small is consequently driven by persistent modelled drug exposure and must not be read as evidence of superior therapy. All parameters are dimensionless and illustrative; no patient data are fitted. Source code and generated tables are included to make every numerical statement reproducible.

Some organisms exhibit declining mortality and increasing fecundity following sexual maturity, a demographic pattern known as negative senescence. According to life history theory, ageing occurs because resources are preferentially allocated to reproduction over somatic maintenance. Models connecting indeterminate growth to negative senescence exist, but none integrate somatic maintenance as a competing allocation decision alongside growth and reproduction. We formulate a life history model in which an individual allocates energy among reproduction, somatic growth, and somatic maintenance and mortality rate depends on both body size and somatic damage. We show that negative actuarial senescence, whereby mortality declines with age, occurs when the proportional change in reproductive value exceeds the proportional change in fitness returns from current investments into reproduction and soma. We derive the necessary conditions for an uninvadable allocation strategy using invasion analysis and Pontryagin's maximum principle, and examine biologically relevant cases numerically. We show that both negative senescence and indeterminate growth arise together as uninvadable outcomes even when maintenance competes for the same resources as growth and reproduction. We show that both diminishing returns to reproduction and diminishing returns to growth can give rise to negative senescence. These results extend the disposable soma theory to organisms with indeterminate growth, in which mortality decreases with size, and identify key mechanisms for the empirically observed association between indeterminate growth and non-senescent demographic trajectories.

Spatial population structure and seasonality are both central to the spread of many infectious diseases of plants, animals and humans. While seasonal forcing in transmission often plays an important role in epidemiological models of a wide range of infectious disease, and we now have some theoretical understanding of the dynamical impacts of spatial structure, the combined effects of these two ubiquitous processes has not been examined in detail. Here, we develop a novel model to explore the combined influence of spatial structure and temporal variability on disease dynamics. Spatial structure is represented using a lattice-based approach with near-neighbour interactions, while temporal variability is included through regular, seasonal, variation of the transmission rate. We use bifurcation analysis of a pair approximation of the full spatial model to identify the parameter regimes associated with qualitatively distinct dynamical behaviours. The model exhibits a remarkably wide range of complex dynamics, including limit cycles, quasi-periodic cycles, multi-year cycles, chaotic dynamics and bistability between these different states. In particular, complex dynamics occur when reproduction is predominantly local, with the dynamics depending critically on the amplitude of the seasonal transmission rate. We show how high transmission rates, high birth rates and in particular low recovery rates are requirements for complex dynamics. We predict that SI-type disease interactions in plant pathogen systems will show complex dynamics even with relatively global transmission dynamics.

The coexistence of competing pathogen strains is shaped by cross-immunity, the cross-protection that infection with one strain confers against another. Although cross-immunity is often asymmetric between strains, this asymmetry is often neglected in the literature on multi-strain coexistence. The effect on coexistence-exclusion outcomes of waning immunity\textemdash which is particularly relevant for antigenically evolving pathogens\textemdash is also poorly understood. To understand how these factors affect strain coexistence, here we analyze a status-based two-strain SIRS model with asymmetric cross-immunity and strain-specific rates for transmission, recovery, and waning of immunity. We derive explicit invasion thresholds that also determine the feasibility and local stability of a unique coexistence equilibrium. Thus, these thresholds allow us to characterize the region of stable strain coexistence, as a function of the cross-immunities and rates of waning immunity. We also obtain closed-form expressions for the strain prevalences at the coexistence equilibrium, showing that the total prevalence may vary non-monotonically as the basic reproduction number of one strain increases. Finally, we show that a transient reduction in transmission can move a coexisting strain pair across an invasion boundary, driving the weaker strain extinct. Applying this result to influenza B, our analysis offers a parsimonious explanation for the disappearance of the Yamagata lineage during the COVID-19 pandemic.

Much of the literature and intuition associated with mathematical epidemiology is driven by deterministic models, which are a reasonable assumption when the population size is large. Stochastic models, especially individual based models, are however considered vital when dealing with small population sizes, especially at times of invasion or extinction. The overwhelming majority of these models (both deterministic and stochastic) assume that the underlying parameters are fixed (or follow a regular seasonal pattern). Here, we consider an analytic framework for dealing with randomly varying parameters through the use of stochastic differential equations - thereby capturing the action of external noisy processes such as weather. In particular, we focus on when the transmission rate, {beta}, varies as the solution to a Cox-Ingersoll-Ross Model, such that {beta} is gamma distributed with autocorrelation. We consider the impact of this parameter variation on a stochastic version of the Susceptible-Infected-Recovered model, and for this 'double-stochastic' model show through simulation and analytical results that exogenous noise increases the impact of stochasticity, potentially leading to more early extinctions, wider variations in the number of cases at equilibrium, but that early growth rate can be faster or slower depending on the precise parameters.

Host-parasite relationships are defined by transmission. Transmission takes diverse forms, with horizontal transmission modes generally characterized as density- or frequency-dependent. However, many observed parasites display mixed-mode transmission using both density- and frequency-dependent routes. To investigate how mixed-mode transmission impacts epidemics, we assessed infection prevalence under single- and mixed-mode transmission when there was a linear trade-off between the probability of infection by the two modes. Our results show that the prevalence of a parasite with mixed-mode transmission can be greater than one with a single-mode only when density- and frequency-dependent routes are associated with different effects of the parasite on host fitness. This work shows that mixed transmission modes per se may not necessarily increase disease prevalence, and that substituting two transmission modes for a single one may result in higher prevalence under limited conditions.

Background: Tuberculosis remains a persistent infectious disease whose control is complicated by latent infection, delayed treatment, incomplete recovery, reinfection, and continuing transmission from infectious individuals. Although treatment is central to tuberculosis management, it is frequently represented only as a transition parameter in mathematical models rather than as a separate epidemiological state. In this study, treatment was therefore incorporated explicitly as an independent compartment so that its influence on transmission, recovery, disease-induced mortality, and long-term disease persistence could be evaluated. Methods: A deterministic nonlinear compartmental model was formulated by dividing the total population into susceptible, exposed, actively infected, treated, and recovered classes. Reinfection of recovered individuals, progression from latent infection to active disease, movement of infectious individuals into treatment, treatment-associated recovery, natural mortality, and disease-induced mortality were included. Positivity and boundedness of the solutions were examined to establish biological validity. The basic reproduction number, R0, was derived through the next-generation matrix approach. Disease-free and endemic equilibria were determined, and their local and conditional global stability properties were investigated using Jacobian analysis, the Routh-Hurwitz criterion, center manifold theory, Lyapunov functions, and LaSalles invariance principle. Normalized sensitivity indices, Latin hypercube sampling, partial rank correlation coefficients, and numerical simulations were also applied. Results: The disease-free equilibrium was shown to be locally asymptotically stable when ,R0<1 whereas sustained transmission and a unique endemic equilibrium were associated with R0>1. Under the stated reduced-model assumptions, stability of the endemic equilibrium was established. Transmission-related parameters were identified as the strongest positive contributors to disease persistence. In contrast, treatment and recovery parameters were found to reduce the reproduction number and infectious burden. Numerical simulations indicated that stronger treatment implementation and reduced transmission opportunities produced substantial reductions in active tuberculosis cases. Conclusion: Treatment was shown to function as both a clinical pathway and an epidemiological control mechanism. The proposed framework may support the design of treatment-centered strategies for reducing tuberculosis prevalence and preventing long-term endemic persistence.

We propose a minimal mathematical framework to describe the go-or-grow dynamics of tumor cells comprising two phenotypically distinct populations. One population is migratory and undergoes constant diffusion, while the other proliferates in an oxygen-dependent manner. The local oxygen concentration governs transitions between these phenotypes. We then ask whether the time scale of cyclic hypoxia, where oxygen is spatially homogeneous but fluctuates between normoxic and hypoxic levels over time, alters the dynamics of the reduced generalist model in the fast phenotype-switching regime compared with the full go-or-grow model. We observe that the reduced generalist phenotype model fails to capture the go-or-grow dynamics under rapid oxygen fluctuations. Our results predict that the tumor growth rate increases in the generalist-phenotype model under rapid oxygen fluctuations compared with tumors of the proposed go-or-grow model. We then explore whether oxygen consumption, which leads to spatial heterogeneity in oxygen concentration, plays a role in distinguishing the dynamics of the reduced model from the proposed go-or-grow model. We demonstrate that the tumor in the proposed go-or-grow model exhibits a higher growth rate than the tumor in the reduced generalist-phenotype model under a slow oxygen consumption rate.

Microbial growth has often been analyzed under the assumption that microorganisms have evolved to optimize phenotypic characteristics of interest, such as growth rate and growth yield. This assumption has been useful, for example, for the genome-scale modeling of metabolism and the study of the allocation of cellular resources to physiological processes. In many experimental situations of interest, however, microorganisms are found to be suboptimal with respect to phenotypic characteristics that are thought to be favorable in that situation. Whereas a strong theoretical framework exists to mathematically relate cellular resource allocation strategies to Pareto optimality of microbial growth and other biological processes, much less is known about the consequences of Pareto suboptimality. We extend the framework to the latter case and show that a given Pareto suboptimal phenotype can be explained by a range of underlying resource allocation strategies, each corresponding to a different growth physiology and biomass composition. We test the predictions with the help of a coarse-grained model of microbial growth and published experimental data, which relate Pareto suboptimal rate-yield phenotypes of Escherichia coli and the microalga Tisochrysis lutea to the macromolecular composition of the cells (storage metabolite and total protein contents). Changing the focus from Pareto optimality to suboptimality provides interesting leads to exploring the diversity of growth strategies that can support a given phenotype. This change of perspective is of practical interest, because some of the growth physiologies within this range may be important for biotechnological applications

A major transition in the evolution of eukaryotic cells was phagocytosis: the ability of one cell to engulf another. In some cases, the engulfed cell persists, forming an endosymbiosis which is often mutually beneficial. However, it is unclear whether the interaction between the ancestors was also cooperative. There are clear benefits to developing the ability to eat a competitor, but can it ever be beneficial to eat a cooperator? We use a dynamical systems framework to identify the ecological conditions under which phagocytosis can evolve from a system of cooperating microbes.Beginning with a system of two cooperating species, we use invasion analysis to identify the parameter regimes under which a phagocytosing mutant can invade and the subsequent population dynamics. We identify a critical threshold that determines whether an invading predator will coexist with the cooperative community or replace members. Repeatedly invading with different predators, we find that the system saturates when there are three coexisting species, though this community does not have a single optimal composition. However, if we incorporate additional costs, we find there is a single evolutionarily stable community of two cooperators and a pure predator. Ultimately, we show that it is possible for phagocytosis to evolve from a cooperative community and, moreover, coexist with the original cooperating species.

Host immunity acquired through prior infections is a key factor shaping the spread of emerging variants. However, accounting for complete infection histories in mathematical models is challenging because the number of possible histories grows exponentially with the number of variants. Here, we present a flexible framework for modeling the co-circulation of multiple variants emerging at arbitrary times while incorporating cross-immunity and variant-specific characteristics. We introduce three tractable approximations: truncation, which retains only infection histories up to a given number of strains; one variant at a time, which is exact when variants circulation does not overlap; and two variants at a time, which is exact when at most two strains co-circulate. We systematically evaluate these approximations across a range of cross-immunity levels and inter-variant emergence delays. Truncation approximations perform best under short emergence delays and strong cross-immunity, whereas the other approaches outperform them at longer delays. When applied to SARS-CoV-2 variant dynamics, the second-order truncation and the two variants at a time approximation accurately capture incidence patterns, whereas the other approaches perform well only for specific subsets of variants. These results identify the parameter regimes in which simplifying assumptions remain valid, provide guidance for model selection, and highlight the key role of infection history in shaping epidemic trajectories.

Persistent viral infections arise from complex interactions between viral replication, host cell responses, and ongoing viral evolution. A general framework linking viral fitness to persistence dynamics is lacking. Here, we develop for the first time a mathematical model of viral persistence that takes into consideration viral fitness variations. Essential to the model is the partition of the classical fitness parameter into three components: replicative, infective, and dispersal fitness. The model was initially inspired by a new experiment on hepatitis C virus (HCV) persistence, established in human hepatoma cells, also reported in this work. This experiment documents two strikingly different viral trajectories depending on the initial replicative fitness of the viral population used to establish persistence. The dynamical model describes the interactions among uninfected cells, infected cells, and infectious virions, and it incorporates, through a continuum, two alternative mechanisms of viral release from cells: budding and lysis. Analysis of the model reveals that viral fitness parameters organise infection outcomes into distinct dynamical regimes. Low replicative and dispersal fitness values lead to viral extinction, whereas high values enable persistence through either stable coexistence or recurrent infection waves. The space of fitness discloses a hierarchy among these parameters, with replicative and dispersal fitness being able to trigger important shifts in the outcome of the infection as opposed to infective fitness. The model identifies trade-offs between replication and dispersal that shape viral production and predicts slow dynamical regimes in which infection may persist despite low detectable viral loads. These dynamical transitions provide candidate mechanisms capable of generating the persistence patterns observed experimentally. Our results establish a computational framework linking multidimensional viral fitness to persistence dynamics and suggest general principles by which evolving RNA viruses transition between extinction (cell curing) and sustained persistence.

Employing DLA methods, this paper explores the self-assembly of collagen fibers and resulting fibrosis at three scales up to the scale of regular lobule models. This allows a mechanistic exploration of the effects of collagen on drug transport (flow and diffusion) and metabolism. In addition, this method permits an analysis of fiber growth characteristics. First, variations of the DLA method of Parkinson et al (1994) will be used to generate multiple explicit collagen microfibril self-assembly using DLA particles in one dimension using cubic grid blocks of (4 mm)3 in a 240 x 20 x 20 grid model. The second stage will be to assess the consequences of various densities of these fibers in three dimensions on flow reductions at a higher scale. Here we utilize DLA methods in cubic grid blocks of (80 nm)3 to mimic 3D collagen self-assembly of fibrils. We then apply a pressure gradient or specified flow rates across a spatially gridded version of these models to quantify flow effects. This region represents a local zone of liver tissue affected by fibrosis. Analytic models of fibrotic effects on flow are employed for comparison. A third stage explores the implications of fibrosis in a liver lobule model using multiple grid blocks of size 3200 mm to represent the lobule tissue. Here, a continuum model of fiber density is employed, based on the previous two scales. The model also includes the effects of additional grid blocks representing sinusoidal flow paths found in the lobule. We contrast and quantify drug propagation and metabolism of molecular dissolved versus nanoparticle delivery vehicles in fibrotic media, achieved by upscaling explicit collagen distributions to appropriate average values.

Buruli ulcer (BU) is a neglected tropical skin disease caused by infection with Mycobacterium ulcerans. The mode of its transmission into the skin has not been fully clarified, but it is believed to involve mechanical transfer, perhaps by insect bites. M. ulcerans secretes a toxin called mycolactone which causes necrosis of the skin and formation of large skin lesions. Understanding the development of these skin lesions over time is crucial to understand the pathogenesis of BU. We formulate a mathematical model to describe the space-time dynamics of M. ulcerans bacteria, mycolactone and skin cell densities as well as fibrin deposition. We investigate the relative contributions to ulcer formation of direct cytotoxicity of mycolactone and ischaemia following the loss of coagulation control. Our results are consistent with existing biological evidence that higher mycolactone concentrations accelerate cell death. Furthermore, the model indicates that the rate of wound expansion is primarily controlled by the diffusion and net birth rate of bacteria, which sets the overall scaling for the speed of wound expansion, along with two dimensionless constants that define the shape of the trajectory of wound edge position against time, namely the ratio of the mycolactone diffusion coefficient to the bacterial diffusion coefficient, and the ratio of the mycolactone degradation rate to the bacterial net birth rate. When mycolactone diffuses significantly faster than bacteria spread or when mycolactone degrades more slowly than bacteria reproduce, the initial development of the wound is faster. In contrast, the long-term wound expansion speed appears to be relatively unchanged across a wide range of the dimensionless parameters.